Glassy Freezing and Long-Range Order in the 3D Random Field XY Model 
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Monte Carlo studies of the 3D random field XY model on simple cubic lattices of size 64 3 , using 
two different isotropic random-field probability distributions of moderate strength, show a glassy 
freezing behavior above T c , and long-range order below T c , consistent with our earlier results for 
weaker and stronger random fields. This model should describe random pinning in vortex lattices 
in type-II superconducting alloys, charge-density wave materials, and decagonal quasicrystals. 



The random-field XY model (RFXYM) was first stud- 
ied forty years ago,[l| as a toy model for random pinning 
effects in the vortex lattice of a type-II superconductor 
in an external magnetic field. For fixed-length classical 
spins the Hamiltonian of the RFXYM is 
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Each fa is a dynamical variable which takes on values be- 
tween and 2tt. The (ij) indicates here a sum over near- 
est neighbors on a simple cubic lattice of size L x L x L. 
We choose each 9i to be an independent identically dis- 
tributed quenched random variable, with the flat proba- 
bility distribution 



P(0i) = 1/27T 



(2) 



for 6i between and 2ir. The probability distribu- 
tion P(hi) is independent of 6^. This Hamiltonian is 
closely related to models of incommensurate charge den- 
sity waves, |2l]3j which include decagonal and other axial 
quasicrystals j3| as special cases. In all of these cases, the 
physical origin of the "random field" is the fact that any 
real sample contains a finite density of point defects, i.e. 
atoms sitting in places where they would not be in the 
ideal structure. 

We will use the term "sample" to denote a Hamiltonian 
with a particular set of the hi and 0i variables. A thermal 
average of some quantity over a single sample will be 
denoted here by angle brackets, ( }, and an average over 
all samples will be denoted by square brackets, [ ]. 

Larkinfl replaced the spin-exchange term of the 
Hamiltonian with a harmonic potential, so that each fa 
is no longer restricted to lie in a compact interval. This 
approximation, which neglects mode coupling, is strictly 
valid only in the limit n — ¥ oo, where n is the number 
of components for each spin. Larkin then argued that, 
for any probability distribution P{hi) which allows non- 
zero values for hi, this model has no long-range ordered 
phase on a lattice whose spatial dimension d is less than 
or equal to four. 

A more intuitive derivation of this result was given 
by Imry and Ma,[f| who assumed that the increase in 
the energy of an L d lattice when the order parameter is 



twisted at a boundary scales as L d ~ 2 . This was later 
justified0,Q by showing that within a perturbative e ex- 
pansion one finds the phenomenon of "dimensional reduc- 
tion" . Term by term for T > T c , the critical exponents of 
any d-dimensional 0(n) random- field model appear to be 
identical to those of an ordinary 0(n) model of dimension 
d-2. 

However, reservations about the correctness of apply- 
ing the n — oo results for finite n, i.e. the substitution 
of the harmonic potential for the exchange term, were 
soon expressed by Grinstein. 0] For the Ising (n = 1) 
case, dimensional reduction was shown rigorously to be 
incorrect, (fil ITol] 

One problem with dimensional reduction for n > 1 is 
that the perturbation theory can only be calculated for 
T > T c , while the twist energy of the order parameter 
is only defined for T < T c . Another, more interesting 
problem for n > 1 is associated with the question of 
properly defining the order parameter. As discussed, for 
example, by Griffiths, jllj there is more than one way to 
do this. For the random- field model, there are several 
different types of two-point correlation functions. Q 

When n = 2, one may use the vector order parameter 
M, defined by 



M 



= N^ 1 cos(fa)x 



sin(fa)y ; 



(3) 



where N = L d is the number of spins. Alternatively, one 
can choose the scalar order parameter M 2 , defined by 
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It is very awkward to calculate M, since its value will 
vary wildly from sample to sample. Thus we expect that 
[M] will be zero, even though (M) for each sample may 
become large at low temperature. On the other hand, we 
anticipate that for weak disorder the probability distri- 
bution for M 2 will become narrow as we take L — > oo. 
Thus it is the quantity (M 2 ), and not the quantity (M), 
which is the natural carrier of long-range order in the 
presence of the random field. 

It is also necessary to consider the possibility that the 
thermal average (M 2 ) for a typical finite sample may be- 
come large as T is lowered, even though (M) remains 
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small. This will occur when the energy required for re- 
orienting (M) is nonzero, but not large compared to T. 
If this were to occur, then there would be no local order 
parameter, and thus no terms in the Landau-Ginsburg 
free-energy functional which depended on gradients of 
(M>. 

Monte Carlo calculations for Eqn. (1) on simple cubic 
lattices were performed some time ago by Gingras and 
Huse0 and FischQ. More recently, additional calcula- 
tions were carried out. [HI The results of these studies, 
which used several different choices for the probability 
distribution P(hi), indicate that for typical samples the 
order parameter (M 2 ) becomes positive at some posi- 
tive temperature T c , as long as the random fields are not 
too strong. However, the precise nature of what hap- 
pens near T c remained unclear. In this work we report 
additional studies using new choices of P(hi), which dis- 
play phenomena which were not seen clearly in the earlier 
studies. 

The computer program which was used was an en- 
hanced version of the one used before. [12J The data were 
obtained from LxLxL simple cubic lattices with L = 64 
using periodic boundary conditions. Some preliminary 
studies for smaller values of L were also done. The 
program approximates 0(2) with Z±2, a 12-state clock 
model. It was modified to enable the study of random- 
field probability distributions of the form 



P(hi) = (1 - p)S(hi - hr) + p5{hi). 



(5) 



This modification makes the program run a few percent 
slower than before. In this work we study two cases. 
Type A samples have h r = 1.5 with p = 0, and Type B 
samples have h r — 2 with p = 0.5. For both of these cases 
we find TJJ near 1.5, which was suggested by Gingras 
and Huse|2| to be optimal for this sort of calculation. 

On a simple cubic lattice of size LxLxL, the magnetic 
structure factor, 5(k), for n = 2 spins is 



L 



" 3 E 



cos(k ■ fij)(cos(<t>i - (f>j)) , (6) 



where ry is the vector on the lattice which starts at site 



i and ends at site j. Thus, setting k to zero yields 
S(0)/L 3 = (M 2 (L)). 



(7) 



As we will see, the extrapolation of (M 2 (L)) to L = oo 
is not trivial. 

In this work, we present results for the average over 
angles of S(k), which we write as S(k). Eight different 
L = 64 samples of the random fields 9i were used for each 
P(hi). The same samples of random fields were used for 
all values of T. 

The Monte Carlo calculations were performed using 
both hot start (i.e. random) initial conditions followed 
by slow cooling, and also cold start {i.e. ferromagnetic) 
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FIG. 1: (color online) Average over 8 samples of ln(S(k)) for 
64 x 64 x 64 lattices with h r = 1.5 and p — at various 
temperatures. The error bars indicate approximately three 
standard deviation statistical errors, and the a;-axis is scaled 
logarithmically. 



initial conditions followed by slow warming. For the Type 
A samples three different cold start initial conditions were 
used, equally spaced around the circle. For the Type B 
samples, where the average random field is somewhat 
weaker, only two cold start initial conditions were used. 

In Fig. 1 we show data for the average over our eight 
Type A samples of ln(5(fc)) for T/J between 1.3125 and 
1.75. The data for T/J > 1.4375 are taken from the hot 
start runs, and the data for lower values of T/J are taken 
from the cold start run which gave the lowest average 
value for the energy, (E). When T — 1.4375 or more, 
all runs for a given sample converge in a rather short 
time to give similar values of (E) and (M) . However, for 
lower values of T/J the hot start runs remain stuck in 
metastable states. 

An attempt was made to drive the hot start runs 
into equilibrium by deep undercooling, followed by slow 
warming. This thermal cycling technique had proven 
successful earlier, 12j for the weak random field case of 
h r = 1.0 with p = 0. For the cases studied here, however, 
a single undercooling was only partially successful in driv- 
ing the samples to equilibrium. The author believes that 
repeating cycling would have worked successfully. How- 
ever, since a detailed study of the nonequilibrium behav- 
ior was not the purpose of the work described here, this 
was not attempted. 

Somewhat surprisingly, the cold start runs were much 
more successful in finding a single free energy minimum 
at T/J between 1.3125 and 1.40625. In the majority 
of samples, all of the cold start initial conditions were 
converging to the same free energy minimum. This was 
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FIG. 2: (color online) Average over 8 samples of ln(S(k)) for 
64 x 64 x 64 lattices with h r = 2 and p — 0.5 at various 
temperatures. The error bars indicate approximately three 
standard deviation statistical errors, and the £-axis is scaled 
logarithmically. 



FIG. 3: (color online) Average over 8 samples of ln(5(7r/64)) 
for 64 x 64 x 64 lattices as a function of T/J. The error 
bars indicate approximately three standard deviation statis- 
tical errors. 



no longer true for T/J < 1.3125, and thus results for 
these lower temperatures are not reported. 

For each sample at each value of T/J, the data shown 
in the figure were obtained by averaging S(k) over 16 spin 
states taken at intervals of 51,200 Monte Carlo steps per 
spin (MCS), after letting the sample equilibrate. Thus 
the length of each data-gathering run was 768,000 MCS. 
The values shown for S(k) are obtained by averaging over 
all of the values of k which have length k. This multi- 
plicity, which is not the same for all values of fc, was not 
taken into account in calculating the error bars shown 
in the figures. Thus each error bar shown in the figures 
represents approximately three standard deviations. The 
procedure used for the Type B samples was essentially 
identical. 

In Fig. 2 we show the corresponding data for the Type 
B samples, for T/J between 1.4375 and 2.00. The data 
for T/J> 1.5625 are taken from the hot start runs, and 
the data for lower values of T/J are taken from the cold 
start run which gave the lower average value for the en- 
ergy, (E). It was also true for these Type B samples 
that the hot start samples were unable to equilibrate on 
accessible time scales for T/J < 1.5625. 

We can define T c in a number of ways. For instance, 
it can be defined as the highest temperature for which 
we have been unable to equilibrate the hot start runs. It 
is also the temperature for which the equilibration time 
for the cold start runs becomes very long. The fact that 
equilibration times get very long at the same T c for both 
the hot start and the cold start initial conditions that the 
behavior of this model can not be explained by the usual 



hysteresis which is associated with a normal first-order 
phase transition. As was discussed in somewhat more 
detail in the earlier work,[l2j the author believes that 
this phase transition should be described by a theory of 
the Anderson- Yuval [l3| type, where the order parameter 
jumps to zero discontinuously at T c . A phase transition of 
this type is also believed to occur in the fc-core percolation 
model, [5] which has been suggested as model for the 
ordinary glass transition. 

The qualitative behavior which we see in Fig. 1 and 
Fig. 2 is that the average over samples of \n(S(k)) shows 
a peak at small k which grows in size as we lower T, 
until it reaches a maximum at T c . The value of T c / J is 
about 1.40625 for the Type A samples and 1.5625 for the 
Type B samples. Upon lowering T below T c , the size of 
the peak becomes somewhat smaller. This is seen more 
clearly in Fig. 3, which shows the temperature depen- 
dence of av[ln(S'(7r/64))] for both types of samples. (The 
value of 7r/64 is approximately 0.049.) 

For T > T c both types display a substantial range 
of T over which S(w/64) is increasing exponentially as 
T is lowered, and this increase slows very close to T c . 
This kind of temperature dependence is well known 
phenomenologically 15] to be a characteristic of glassy 
freezing behavior. However, the reader should keep in 
mind that one expects the usual glass transition to cor- 
respond to n = 3 spins, and not the n = 2 spins studied 
here. If one generalizes Eqn. (1) to the n — 3 case, sub- 
stantially more computing resources would be required 
to make the equivalent calculation feasible. 

Below T c , 5(7r/64) first decreases, and then becomes 
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approximately constant, at least over the range of T 
shown in the figure, where the samples are believed to be 
equilibrated. This behavior below T c is consistent with 
the existence of a 5-function peak in S at k — in the 
limit L — > oo, indicative of true long-range order. One 
would not expect to see this behavior if there were only 
power-law correlations in S"(fc), with no <5-function peak. 

The reader should understand that we are not claim- 
ing the existence of some mathematical error in the 
analysis [3| of the elastic glass, which leads to the power- 
law correlated Bragg glass phase. That analysis begins 
with the Larkin Hamiltonian, [lj which approximates the 
spin-exchange term of Eqn. (1) with a harmonic interac- 
tion. What we are saying is that the long-range behav- 
ior of the Larkin Hamiltonian is different from that of 
Eqn. (1), at least for the types of P(hi) distributions we 
have studied. 

If we fit the data for av[ln(S'(fc))] at T c over the range 
7r/64 < k < 7r/8, as shown in Fig. 1 and Fig. 2, we 
find good straight line fits, with slopes of -3.20(3) for the 
Type A samples and -3.10(2) for the Type B samples. 
Because both of these slopes are more negative than -3, 
it is clearly not correct to extrapolate these fits to k = 0. 
S(k) satisfies a sum rule, and thus nonintegrable singu- 
larities cannot occur. Thus if we were able to obtain 
data using these P(hi) distributions for even larger val- 
ues of L, we would expect to find that the behavior of 
av[ln(S'(fc))] at T c for small enough nonzero k would look 
like the data for the case h r — 2 with p = found in the 
earlier work, fljj] 

In this work we have performed Monte Carlo studies 
of the 3D RFXYM on L = 64 simple cubic lattices, with 
isotropic random field distributions of two types: h r = 
1.5 with p = 0, and h r = 2 with p = 0.5. We present re- 
sults for the structure factor, S(k), at a sequence of tem- 



peratures. In agreement with prior results 12] on weaker 
and stronger random fields, we argue that our results in- 
dicate a phase transition into a long-range ordered low 
temperature phase, with a jump in the order parameter 
at T c . The behavior above T c strongly resembles well 
known phenomenology for glassy phase transitions. 15 
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